function [L D] = ldl_chol(A)
    if norm(A - A') > eps
        fprintf("The input matrix must be symetric.\n");
        return;
    end
    n = size(A, 1);
    d = zeros(n, 1);
    L = zeros(n);
    C = zeros(n);
    for j = 1:n
        dl2 = 0;
        for s = 1:j - 1
            dl2 = dl2 + d(s) * L(j, s)^2;  
        end
        C(j, j) = A(j, j) - dl2;
        d(j) = C(j, j);
        for i = j + 1:n
            dll = 0;
            for s = 1:j - 1
                dll = dll + d(s) * L(i, s) * L(j, s);
            end
            C(i, j) = A(i, j) - dll;
            L(i, j) = C(i, j) / d(j);
        end
        L(j, j) = 1;
    end
    D = diag(d);
end